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Abstract 

We report on progress in algorithms for iterative phase retrieval. The theory of convex optimization 
is used to develop and to gain insight into counterparts for the nonconvex problem of phase retrieval. We 
propose a relaxation of averaged alternating reflectors and determine the fixed point set of the related 
operator in the convex case. A numerical study supports our theoretical observations and demonstrates 
the effectiveness of the algorithm compared to the current state of the art. 



1 Introduction 

The phase retrieval problem is a classical inverse problem in optics that has received renewed interest in 
applications to nonperiodic scatterers and macromolecules. While scattering from some structures allows 
one explicitly to compute the phase from magnitude measurements [12], more general classes of scatterers 
require the use of less direct methods. So called iterative transform methods pioneered by Gerchberg and 
Saxton [11], and Fienup [10] are well established generic techniques for iteratively recovering the phase in a 
variety of settings. Recent developments in imaging [9, 13, 14, 18-21], have placed a premium on improving 
the efBciency and stability of phase retrieval algorithms. 

In this work we derive a stable and fast new strategy for phase retrieval, what we call Relaxed Averaged 
Alternating Reflection (RAAR), that falls under the category of iterative transform methods [15]. The 
motivation for the RAAR algorithm comes from recent work in which another new algorithm, the Hybrid 
Projection Reflection algorithm (HPR), was presented [4]. The HPR algorithm was originally concieved as 
a single parameter relaxation of the well known Douglas-Rachford algorithm applied to phase retrieval. The 
HPR algorithm can also be viewed a special case of the three parameter difference map recently proposed 
by Elser [7]. 

There are two fundamental and distinct issues that accompany these algorithms. The first is the incor- 
poration of a priori information into the constraint structure of the algorithms. The second is the choice 
of algorithm parameter values. Regarding the first issue, it is difficult to overestimate the effect of the 
constraints on the mathematical properties and performance of the algorithms. There have been several 
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studies on the choice of constraints in apphcations to crystaUography [8, 16, 17]. We use a simple example 
to illustrate how seemingly minor changes in the physical domain constraints can lead to algorithms that 
appear very different. This has caused some confusion in the literature which we hope to clarify through 
an examination of the abstract algorithmic structures behind the leading techniques. The choice of param- 
eters also has a dramatic impact on the mathematical properties of the algorithms and hence performance. 
Physical insight often provides the best (and only) basis for chosing values for the algorithm parameters, but 
this is not always available or reliable. In the case of the HPR algorithm, our numerical experiments have 
not provided an empirical basis upon which to make recommendations. A more mathematically rigorous 
approach also appears to be difficult and has been found in only a few very special cases. For instance, in the 
convex setting the convergence properties of the HPR algorithm are known for the unrelaxed case [5] . For 
the relaxed HPR algorithm and the more general difference map a complete and mathematically rigorous 
analysis has yet to be found. To circumvent the analytical barriers facing the difference map and the HPR 
algorithm, we introduce the RAAR algorithm, which is conceptually simple, analytically tractable and easy 
to implement; moreoever, it outperforms the current state of the art. While the RAAR algorithm coincides 
with the HPR algorithm in a limiting case, it does not fall in the class of algorithms covered by Elser's 
difference map framework. 

A precise statement of the leading algorithms is given in Section [21 In this same section we provide a 
terse outline of the mathematical justification for the RAAR algorithm. In Section O we demonstrate the 
effectiveness of the algorithm and make practical recommendations for implementation. 



2 Phase Retrieval and Iterative Transform Algorithms 
2.1 Phase retrieval 

We are interested in recovering the scattering amplitude of a medium that has been illuminated by an 
electromagnetic wave from measurements of its spatial coherence function and other a priori information. 
For the sake of concreteness, we assume that is a real-valued, nonnegative function supported on some 
prescribed bounded set D, that is £ 9 m* : ^ R+ with supp(u*) C C Z^. Here C is a Hilbert 
space of square integrable functions, is the domain - in this case the physical domain - corresponding 
to discrete (i.e. sampled) waves, R+ is the positive reals and supp(u,) is the support of u,. Writing this in 
terms of constraints, we have u* G 5+ C where 5*+ is the set of nonnegative functions in C with support 
on D. If we require only that the functions be supported on £), we denote the corresponding constraint 
set by S. The sets S and 5*+ are refered to as the physical domain constraints. The other constraint we 
consider comes from the the data, m, which we presume consists of noisy magnitude measurements in the 
far field, thus m is proportional to the modulus of the Fourier transform of it,. We therefore refer to the 
domain of the image data m as the Fourier domain. In terms of constraint sets, we write that G M where 
M = {w G £ I |jFw| — m} and J-v denotes the discrete Fourier transform of v. We shall refer to the set M 
as the Fourier, or image domain constraint. Note that 5+ is a convex set, while M is nonconvex. It is the 
nonconvexity of the magnitude constraint that does not allow us to transfer classical convergence results for 
the most common algorithms to the case of phase retrieval. For further discussion see [3]. 
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2.2 Iterative Transform Algorithms 



We formulate the problem of phase retrieval as a feasibility problem: 

find ue S+n M. 

Iterative transform techniques are built upon combining projections onto the sets 5+ and M in some fashion. 
While they are seldom written as fixed-point algorithms, iterative transform algorithms can usually be put 
into the form u„+i = Tu„ where T is a generic operator in which the projections and averaging operations are 
embedded (see [3,4]). For added control and flexibility, one often includes a relaxation strategy parameterized 
by (3. We write the relaxed operator with generic, single parameter relaxation strategy V (there can be 
infinitely many such strategies) as V(T, /3). In order effectively to exploit relaxations for improved algorithm 
performance it is necessary to understand the mathematical properties of the operator V(T, /?) - first and 
foremost of these is the characterization of the set of fixed points. Fix V(T, /3). We return to this issue at 
the end of this section. 

The operators we study are built upon projectors and reflectors. Denote by Pc an arbitrary but fixed 
selection, or projector, from the possibly multi-valued projection onto a subset C of C. Closely related is the 
corresponding reflector with respect to C 

where / is the identity operator. By definition, for every u ^ C, Pc{u) is the midpoint between u and Bc{u)- 
Specializing to our application, the projector, PmU, of a signal u ^ C onto the Fourier magnitude constraint 
set M is given by 

Pm{u)^T \vo), where wo(0 = < K"(OI (1) 

[m(^), otherwise . 

Here, J-^^ is the discrete inverse Fourier transform and a- selection from the multi- valued Fourier domain 
projection. For further discussion of this projector see Luke et al [15, Corollary 4.3] and [6]. We return to 
the issue of niultivaluedness of the magnitude projection in Section |3| The projection of a signal u G £ onto 
5"+ is single- valued (since S+ is convex), and is given by 

(V.eZ^) (P,,(^.))(.) = h""^°'"("^^' ll^^2e 



One of the best known iterative transform algorithms is Fienup's Hybrid Input-Output algorithm (HIO) 
[10]. We use this as our benchmark for performance. In the present setting, HIO is given as 

^ (\ \{PM{un)){x), \i X e D &nd {P^,{un))(x)>{)] 

\un[x) - Pn[PM{Un)j{x), Otherwise. 



There have been several attempts to identify the HIO algorithm with a broader class of relaxation strategies 
that can be written as fixed point iterations, that is, in the form Un+i = V{T , [3n)un.. Bauschke, Combettes 
and Luke [3] proved that, when only a support constraint as opposed to support and nonegativity is applied 
in the physical domain, then the HIO algorithm with jS — \ corresponds to the classical Douglas-Rachford 
algorithm for which convergence results in the convex setting are well known. In a subsequent article 
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Bauschke Combettes and Luke [4] proved that, for physical domain support constraints only, the HIO 
algorithm corresponds to a particular relaxation of the Douglas-Rachford algorithm, that is 

[u„(xj - Pn[-rM[Un))[x), othcrwisc, 

is equivalent to 

Un+l = ^{Rs{Rm + iPn - 1)Pm) +/+(!- /9„)Pm) (5) 

Independent of these results, Elser [7] showed the correspondence between the HIO algorithm with only 
support constraints in the physical domain and the difference map, 

Un+l = (I + P {Ps ((1 - 12)Pm - 72/) + Pm ((1 - ll)Ps - 111))) K), (6) 

for the case where 71 = —1 and 72 = The correspondence between the difference map and the HIO 

algorithm does not carry over to the case of support and nonnegativity constraints. The correct formulation 
of the corresponding algorithm was given in [4, Proposition 2], where it is shown that 

Un+l - \ {Rs^ (Rm + iPn - 1)Pm) + / + (1 - /3„)-Pm) (7) 

is equivalent to 

{{PM{un)){x), if X e and 

(i?A/K))(a;) > (1 - l3n){PM{un))ix); (8) 
Un{x) ~ Pn{PM{un)){x), Otherwise. 

In [4] the fixed point iteration {Tj) is called the Hybrid Projection Reflection (HPR) algorithm, which is 
equivalent to the difference map (with 71 = — 1 and 72 = 1//3) applied to support and nonnegativity 
constraints: 

Un+l + {Ps+ ((1 - 72)Pm - 72/) + Pm ((1 - ll)Ps^ - 111))) (U„). (9) 



It is important to note that, while the form of prescriptions of projection algorithms in terms of fixed 
point iterations Un+i — V(T, Pn)un does not depend on the underlying constraints, this is not the case for 
prescriptions of the form © , Q and As we have seen, slight changes in the constraint sets can result in 
dramatic changes in the form of algorithms when written in this way. When written as fixed point iterations, 
the effect of changing the constraint structure is seen in the mathematical properties of the operator rather 
than the form of the algorithm. 

Preliminary numerical results indicate that the HPR algorithm is a promising alternative to HIO - HPR 
is more stable and, at least with simulated data, produces higher quality images. Detailed convergence 
results have been obtained in [5] for the unrelaxed HPR algorithm (/3 = 1) in a convex setting. At this time, 
however, there are no mathematically rigorous results proving convergence or suggesting how to choose the 
relaxation parameter /3„ to improve performance. Another drawback to the HPR algorithm is that, while it 
consistently delivers higher quality solutions than HIO, it can take longer to achieve this. The algorithm we 
propose next addresses both the analytical drawbacks as well as the performance issues regarding the HPR 
algorithm and the more general difference map. 

The new algorithm we propose is given by the following: given any uq G £ , generate the sequence 
wo,ui,U2, ... by 

Un+l=V{%,l3n)Un (10) 
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where 

V{%,(3)^P% + {\-P)P,, and % = \{Rs^Rm+I). (H) 

To underscore the connection of this algorithm with the Averaged Alternating Reflection (AAR) algorithm 
studied in [5] , we refer to IjlUI) as the relaxed averaged alternating reflection (RAAR) algorithm. For [3 — 1 
the RAAR, HPR, AAR, and the difference map (71 = —1 and 72 = 1//3) algorithms are equivalent. For 
/3 7^ 1 the RAAR algorithm is fundamentally different than HPR; moreover, it cannot be derived as a special 
case of the difference map The recursion (|10|) can be written analogously to ^ and To sec this, we 
proceed as in Proposition 2 of [4]. Given an arbitrary signal v ^ C, let w"*" = max{w,0} and — min{u,0} 
be its positive and negative parts, respectively. Then can be rewritten as 

Un+i = [-Xd^ ■ /3„(2Pm - /) - [Xd ■ /3„(2Pm - I)V + Pm) M- (12) 

There are 3 cases to consider: (i) If x G D and {RMUn){x) > 0, then ifT^ yields Un+i = Pm] (h) ii x ^ D and 
{RMUn){x) < 0, then becomes u„+i(a;) = ( ((1 - 2f3n)PM + Pnl) {un)){x)] (iii) if t ^ D, then can 
also be written as u„+i(a:) = ( ((1 — 2/5„)i^f + /?„/) {un)){x). Altogether this yields the following algorithm 

f iPM{un)){x), iix(zD and (R,,{un)){x) > 0; 

(VxGZ^) u„+i(x) = r (13) 

[PnUnix) - (1 - 2f3n){PM{un))ix), othcrwisc. 

We summarize the above discussion in the following proposition. 

Proposition 2.1. Algorithm is equivalent to the recursion ((TUIl. 

The update rule in algorithm ((T^ depends on the pointwise sign of the reflector (i?M(un)) (a^) whereas 
the update rule for Fienup's HIO algorithm ^ depends on the pointwise sign of the projector {PM{un)){x). 
The difference between the RAAR update rule and that for HPR © is much starker. Also note that the 
"otherwise" action is simply a relaxation of the conditional action in the HIO algorithm; this is, again, very 
different than the HPR algorithm. 



2.3 The RAAR algorithm: convex analysis 

To gain some insight into the behavior of the algorithm above, we study the behavior of the convex analog to 
V{%,(3). Let A and B be two closed convex subsets of C. Replace S+ and M hy A and B respectively. Let 
E C A denote the set of points in A nearest to B, and let F C i? denote the set of points in B nearest to A. 
The gap vector between A and B, denoted hy g G £, is defined by 5 = Pc\{B-A}(^y Loosely interpreted, this 
is a vector pointing from E to F with \\g\\ measuring the smallest distance between A and B. For instance, 
it An B then g = 0. For a more precise treatment see [1,2]. The convex counterpart to (|ll|l . the central 
operator in the RAAR algorithm, is defined by 

V{T,,(3) = f3T, + (1 - (3)Ps, < /3 < 1 where = ^{R^Rs + T). (14) 

When discussing convergence of projection-type algorithms, one must take care to distinguish between 
consistent and inconsistent feasibility problems. In the current convex setting, consistent problems satisfy 
A n B ^ 0; when An B — the problem is said to be inconsistent. Inconsistent problems are common in 
applications where the a priori information represented by the constraint sets is highly idealized, particularly 
in the presence of noise. Bauschke, Combettes and Luke [5] show that the properties of the AAR algorithm 
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(that is, RAAR with /3 = 1) for consistent problems are very different from inconsistent problems. The 
reason for this is that the operator T* does not have a fixed point ii An B — 0. For < /? < 1 the convex 
instance of the RAAR algorithm avoids these complications by transferring questions of consistency of the 
constraints to the existence of nearest points. In other words, the RAAR operator enjoys the advantage that 
Fix F(r»,/3) is independent of whether or not the associated feasibility problem is consistent. This is the 
content of the following theorem. 



Theorem 2.2. LetO < f] <l. Then 



FixV{n,/3)^F~^g (15) 



where g is the gap vector between A and B and F C B is the set of points in B nearest to A. Moreover, for 
every u S Fix V{T^, (3), we have the following: 

(i) u = PbU -g\ [ii) PgU — P^RgU = g; {Hi) PgU G F and PaPbU G E. (16) 

1 fi 

Proof. To prove the result we must show (a) that F — /3g/{l ~ /3) C Fix V{T^,, (3) and (b), conversely, that 

Fix V^(T»,/3) C F ~ I3g/{1 — /?). The first statement (a) is proved analogously to the proof of equation (18) 
of [5] . In the interest of brevity, we leave this as an exercise. 

We show that Fix(/3T, + (1 - /3)Pb) C F ~ j^g. To see this, pick any u e Fix(/3r, + (1 - P)Pb). Let 
f — PbU and y = u — f. For any b E B, since i? is a nonempty closed convex set and f — PbU, we have 
{b — PbU, u — f) < 0. which yields 

{b-f,y)^{b-f,u-f)< 0. (17) 
Recall that iA(2/ — u) — Pa{2PbU — u) — PaRbU. Together with the identity [5, Proposition 3.3(i)] 

(Vu G £) u~T^u = PbU - PaRbU (18) 

equation H17|l yields 

P^(2/-u) = / + T,u-u. (19) 

Now l3T^u + (1 - I3)PbU = u yields 

T^u~u^^^-—^{u~Pbu). (20) 

Then (jTHJl and give Pa(2/ - u) = / + ^-^{u - f) = f + ^-j^V- As above, for any a e A, since A 
nonempty, closed and convex, we have (a — J^(2/ — u), (2/ — u) — ^^(2/ — u)) < 0, and hence 

> (.-(/ + i^.),(2/-.)-(./+i^. 

= ( « - ./ H — 3— y h -y 3— y 



IS 



l(-„ + ;.,) + <i_«W^ (21) 



Now ((T7|) and yield {b — a, y) < — < 0. Now take a sequence ap, ai, 02, ... in A and a sequence 





bo,bi,b2, ■ ■ ■ i^o. B such that g„ — 6„ — a„ g. Then 



(VneN) {gn,y)< ^||yr<0. (22) 
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Taking the limit and using the Cauchy-Schwarz inequality yields 



|y||< Y^llffl!- (23) 



Conversely, u - {fiT^u + (1 - (3)Psu) ^ (3 {f - P^(2/ - u)) + (1 - /3)j/ = gives 



\y\ 



1-/3 



f-PA2f-u) >Y^ll5l|. (24) 



Hence ||y|| = y^||w|| and, taking the limit in H22() . y = —j^g, which confirms (i). It follows immediately 
that / — PaRbU = g which proves (ii) and, by definition, implies that PgU ^ f Cz F and P^PbU £ E. This 
yields (iii) and proves H15|l . ■ 

In words, regardless of whether or not AoB is empty, as long as there are points in B that are nearest to A, 
then the RAAR operator V{T^, (3) has a set of fixed points, and these are precisely the points in B nearest 
to A, translated by the scaled gap vector. This is the starting point for the convex heuristics behind the 
RAAR algorithm. Statements about convergence and more detailed behavior of the algorithm are beyond 
the scope of this work. 

We conclude the mathematical analysis with some observations that motivate the relaxation strategy we 
implement in Sectional We wish to use the parameter (3 to control the step size between successive iterates 
and, as much as possible, to steer the iterates. Far away from the solution, it is easy to see the damping 
effect of the parameter < /? < 1, which derives from the form of the relaxation (|11() as simply a convex 
combination of the operator T^, and the projector onto the data Pm - the smaller the relaxation parameter /3, 
the closer to the data we require the iterates to stay. It was noted in [4] that, regardless of the relaxation, the 
HPR algorithm (jHl takes significantly longer than the HIO algorithm © to reach a suitable neighborhood of 
the solution, although, once near a solution, HPR delivers consistently better images with greater stability 
and reliability than HIO. We show in the next section that the dampening effect of the relaxation in the 
RAAR algorithm is just what is needed to control the initial behavior of the HPR algorithm. 

For the behavior of the algorithm near the solution, we rely on the convex analysis. By (|15|l . the relaxation 
parameter /3 effects the fixed points of the operator through the gap vector. If the feasibility problem is 
consistent, that is, A n -B ^ 0, then the gap vector 5 = 0. In this case, is it not clear what effect, if any, (3 
will have on convergence. On the other hand, if the problem is inconsistent, that is, A n B = 0, and 5 0, 
then, by 1)15(1 , the set of nearest points F can be translated arbitrarily far away in the direction g by letting 
(3 approach 1 from below. We use this to gain some control on the step size between successive iterates and 
the directions of the steps. 



Proposition 2.3. Let m„ G C satisfy ||u„ - < 5 where up^ e Fix V"(T,,/3„) and < /3„ < 1. Define 
Un+i = V{T^, Pn+i)un for any < /3„+i < 1. Then 



Un+l - fp. 



'n+1 



1-A 



-g 



< 6, where fp^ = PsUp^^ G F. 



(25) 



Proof. For any m G £, by lO, we have T^(r*, /3„+i)u - F(T*,/?„)u = (/3n+i - Pn) {Pa - I) RbU, which, 
together with ifT^ fil. yields 



l-f3n 



l-f3n' 



Since V{T^, Pn+i) is nonexpansive, the result follows from (|26|l . 
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While the HPR algorithm gives quite stable solutions eventually, the above theory suggests that this stability 
can be improved in a controlled fashion. Consider the fixed point iteration as a descent algorithm minimizing 
some error metric (in fact, minimizing the gap distance) where —g is the direction of descent. By (|25(l and 
the first equation in H26|) . 

Un+i W V{I^,Pn+l)Uf3„ = — g, 

t Pn 

thus one can use /3„+i to affect steps in the direction —g ranging, in the limit, from length — /?„/ (1 — /3n) to 
1 as (3n+i varies from to 1 respectively. The difference m„+i — u„ for the unrelaxed algorithm {(3 = 1) was 
shown in [5] to converge to the negative gap vector — g in the inconsistent case. The effect of the relaxation 
is primarily to dampen the iteration in the neighborhood of a solution in the case of inconsistent problems. 
To see the advantage of this, consider the nonconvex case and suppose that the problem is inconsistent (that 
is, the gap vector g 0). The only case of the HPR algorithm for which we can say anything is the case 
(3=1, which is the same as the unrelaxed RAAR (or AAR) algorithm, so we restrict the discussion to the 
RAAR and AAR algorithms. The convex analysis of the AAR algorithm shows that, even though the gap is 
attained, the iterates m„ continue to move in the direction —g without end. In the nonconvex setting, even 
if the true gap is attained, the continued progress of the iterates in the direction —g could push the iterates 
away from the domain of attraction of the local solution and into a different domain of attraction. Thus the 
projections of the iterates, or the shadows might never converge. This "wandering" of the iterates near an 
apparent local solution has been observed both with the HIO and HPR algorithms, though it is much less 
severe and destabilizing with HPR than it is with HIO. The relaxations in the RAAR algorithm can be used 
to either dampen the iterates near a local solution to slow drifting out of a domain of attraction, or to halt 
the wandering of the iterates altogether by holding the relaxation parameter at a fixed value less than 1. 



3 Numerical Implementation 



Our goal with the RAAR algorithm is to use dynamic relaxations to shorten the initial "warm-up" phase of 
the HPR algorithm and to stabilize the algorithm near a local solution. The algorithm we consider is 

Un+l ~V{%,l3n)Un. (27) 



Before outlining our specific implementation, some remarks are in order about the calculation of T given 
by As discussed in [15, Section 5.2] the projection onto the magnitude constraint is a numerically 

unstable operation due to the multivaluedness of the projection operator. We therefore recommend the 
following approximation to Pm (see [15, Eq.74]): 

P.,u . VJ.u - m] ^^^^^^Tu] (28) 

VVd-^Wp+e')^ ; £2)3/2 J 



for < e ^ 1 , where 



1 /„ „9 ll2\ , _ ITul'^ 



Mu)= -(\\u\\^-\\T-'d-m\n, where V = ^ (29) 

(l^«l'+62) 

Define 

V{%, /3) = i {R,^ (2V J, -/)+/). (30) 
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Figure 1: Original images and corresponding data used for the comparison of the HIO and IfPR algorithms. The 
center of (a) is a 38 x 38 pixel section of the standard Lena image, zero-padded to a 128 x 128 matrix. Frame (b) 
is the noiseless Fourier magnitude data m corresponding to image (a). The same object domain support constraint 
(and initial guess) of size 64 x 64 pixels, shown in (c), is used for each trial. 




Under reasonable assumptions, by the continuity of Rs^ and [15, Corollary 5.3] it can be shown that 
VJe(M)^P„(w) and V{%,f])u-^V{%,f])u as e 0. 

Using the stable approximation V{%,, {3) given by (|30|l . from the initial guess uq we generate the sequence 
Mo,Mi,M2, ■ • ■ by 

zi„+i =l/(7;,/3„)un where /3„+i = /3o + (1 - /3o) (l - exp (-(n/?)^)) . (31) 

The rule for updating /3„ is a smooth approximation to a step function from the value /?o to the value 1 
centered at iteration n — 7. We compare this algorithm to the HIO (O and HPR ((SJ algorithms using the 
same stable projection approximation. We study algorithm performance with noisy data. The initial points 
Mo are chosen to be the normalized characteristic function of the support constraint shown in Figure ^^c). 

The data consists of the support/nonnegativity constraint, shown in Figure^c), and Fourier magnitude 
data m, shown in Figure ^b), with additive noise rj - a symmetric, randomly generated array with a zero 
mean Gaussian distribution. The signal-to-noise ratio (SNR) is 201ogiQ ||m||/||M|| — 34 dB. As motivated 
in [4], the error metric we use to monitor the algorithms, £b^, is given by 

Es^ {Xn ) = " , " „2 ^" . (32) 

||-Pm(w„)|| 

We compute the mean value of the error measure Es+ over 100 trials with different realizations of the noise 
and the same initial guess. 

First, we compare the mean behavior over 100 iterations of two sets of realizations of the algorithms, each 
corresponding to different relaxation strategies, (3 — 0.75, (i — 0.87, (3 = 0.99 and variable /3„ governed by 
()31|) with — 0.75. The average value of the error metric at iteration n, Es^{xn), is shown in Figure [21 
These are all given in decibels (recall that the decibel value of a > is lOlogig(a)). In Figure 13 we show 
typical estimates generated by the respective algorithms at iteration 35, all from the same realization of noise 
and the same initial guess. While the RAAR algorithm with /? = 0.75 appears to perform well as measured 
by £^3^ (see FigureEIa)), it is clear from Figure|2|that the quality of solutions found by the RAAR algorithm 
degrades rapidly as the relaxation parameter /3 becomes small. For values of f3 near 1.0 the quality of the 
iterates generated by the RAAR algorithm does eventually improve, however, as with the HPR algorithm, 
it takes many more iterations to achieve this imporvement. For static values of P the best performance for 
the RAAR algorithm appears to be achieved with a value of /? = 0.87. The variable /?„ trials for the RAAR 
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Figure 2: Error metric Es^{x„) averaged over 100 realizations of noise (SNR=34 dB). For (a)-(c) the relaxation 
parameter for the respective algorithms, /3„, is fixed. For (d) /3„ varies from 0.75 to 1.0 according to 13111 . 



P =0.75 p =0.87 




algorithm yielded the best overall results, measured both by the error metric, as well as observed picture 
quality. In contrast to this, the relaxation parameter does not appear to have any identifiable effect on the 
performance of the HIO or HPR algorithms. 



4 Concluding Remarks 

There are infinitely many relaxation strategies one could implement for iterative transform methods, but 
very few of them admit a meaningful mathematical analysis. The standard for phase retrieval algorithms, 
Fienup's HIO algorithm, has been identified in a special case with the promising HPR algorithm, which in 
turn, has been identified as a special case of Elser's difference map. For each of these algorithmic frameworks, 
the mathematical properties of the algorithms vary drammatically with the parameter values in a manner 
analogous to bifucations of dynamical systems. A complete mathematical analysis must treat all relevant 
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Figure 3: Typical images recovered after 35 iterations of the HIO, HPR, and RAAR algorithms for different 
relaxation strategies with the same realization of data noise (SNR=34 dB) and the same normalized initial guess. 
The variable /3n trials were generated according to the rule given by 13H . 



HIO HPR RAAR 
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intervals of parameter values on a case by case basis. No such analysis is available for the HIO, HPR or 
difference map algorithms. To circumvent these difficulties and to improve upon the HPR algorithm, we 
propose a simple relaxation, the RAAR algorithm, of a well understood Averaged Averaged Reflection (AAR) 
algorithm. The relaxation is a convex combination of the AAR fixed point operator, and the projection onto 
the data. This intuitive framework is mathematically tractable and provides an easy strategy for the choice 
of relaxation parameter that, moreover, improves algorithm performance. In contrast, it appears that similar 
relaxation strategies have little effect on either the HIO or the HPR algorithm. We cannot suggest a rule 
by which to select a static value oi (3 - this depends on the data. Nevertheless, based on the results for 
the variable /3„ trials, we can recommend the fairly generic dynamic relaxation strategy of 1)31(1 for getting 
the best performance from the RAAR algorithm. Here the algorithm is significantly relaxed in the early 
iterations, helping the algorithm quickly to find a neighborhood of the solution while maintaining fidelity 
to the data, and then decreasing the relaxation (i.e. increasing /3„) in the neighborhood of the solution to 
avoid stagnation at a poor local minimum. To stabilize iterates in the domain of attraction of a solution, a 
final fixed value of /3 close to, but less than, 1, say P — .99999 should be chosen. In a technical point, we 
also proposed a smooth perturbation of the magnitude projector H28|) to improve the numerical stability of 
computing the projection onto magnitude constraints. 
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